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Summary 


A numerical technique to solve the Euler equations for steady, 
one-dimensional flows is presented. The technique is essentially 
implicit, but is structured as a sequence of explicit solutions 
for each Riemann variable separately. Each solution is obtained 
by integrating in the direction prescribed by the propagation of 
the Riemann variables. The technique is second-order accurate. 
It requires very few steps for convergence, and each step re- 
quires a minimal number of operations. Therefore, it is three 
orders of magnitude more efficient than a standard time-dependent 
technique. The technique works well for transonic flows and pro- 
vides shock fitting with errors as small as 0.001. Results are 
presented for subsonic and transonic problems. Errors are 
evaluated by comparison with exact solutions. 


_1_. Introduction 

Numerical solutions of the Euler equations, including accu- 
rate calculations of discontinuities, should provide the best 
description of compressible, inviscid flows. Nevertheless, when 
such flows are steady, a strong emphasis has been put, in recent 
years, on solvers of the potential equation, which apply relaxa- 
tion methods. The assumption of a velocity potential inhibits 
variations of entropy. Consequently, positive jumps in pressure, 
commonly called "shocks", occur at constant entropy and do not 
satisfy the Rankine-Hugoniot conditions; legitimate doubts arise 
on the accuracy of the results. Strong arguments, however, have 
been set forth on defense of potential-relaxation techniques; 
above all, their speed of execution, and the possibility of ex- 
ploiting multiple grid or multi-grid accelerating effects. Wha- 
tever accuracy may be obtained when vorticity and entropy effects 
are neglected, it is aimed at by using a very large number of 
computational nodes. 

We claim that a good Euler solver can reach an equal (or 
higher) degree of accuracy with much fewer nodes. Recently, sub- 
stantial progress has been made in the numerical integration of 
the Euler equations, by introducing robust and accurate integra- 
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tion schemes [l, 2,3, 4, 5], methods for the correct handling of 
boundary conditions and means for the definition of suitable com- 
putational grids (these topics are reviewed in [6] and [ 7 ], 
respectively). Unfortunately, if a steady state solution has to 
be reached using a standard time- dependent Euler solver, conver- 
gence to the asymptotic solution is slow, and the number of 
operations to be performed per step is larger than in a 
potential-relaxation procedure. Consequently, the calculation 
turns out to be much more expensive, despite the use of a coarser 
grid. 


We present here an extremely simple Euler solver for steady 
flows, which requires a minimal amount of operations per step and 
a small number of steps to achieve convergence. In principle, it 
amounts to solving two bidiagonal matrices per step. 

For a better illustration of the method, we limit the 
present report to one-dimensional problems. We will examine our 
findings on two-dimensional problems in a successive paper. 

The technique presented here is derived directly from the 
same ideas which produced the "lambda ,, -scheme in its more recent 
formulation [8]. Conscious and unconscious borrowing of ideas, 
however, is part of the creative labor of every new project and, 
whenever possible, should be acknowledged. The idea of sweeping 
(that is, of attempting to solve a bidiagonal matrix instead of a 
tridiagonal one) was inspired by MacCormack [ 9 ]. The modeling of 
boundary conditions is essentially due to M. Pandolfi, and L. 
Zannetti should be thanked here for his active faith in Riemann 
variables and an early suggestion of the relation between AR^ a nd 
the shock Mach number. I wish also to acknowledge the collabora- 
tion of Miss Catherine Fahy in the preparation of the code and 
the organization and analysis of data, on a level above clerical 
routine. 


2. Outline of the technique 

To explain the origin of the technique, let us begin by wri- 
ting the Euler equations in characteristic form: 
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H 1t - as t + ^ 1 (R 1 x - as*) + d « 0 


R 2t " aB t + X 2^ R 2 x “ a8 x > + d ” 0 (1) 

8 t + u 3 x " 0 

where 

R 1 * a/S - u, R 2 “ a/5 + u (2) 

X^ - u - a, Xg - u + a (3) 

d - a u A x /A (4) 

and a, u, a are the speed of sound, particle velocity and entro- 
py, respectively; 5 » (y-1 )/2, and y is the ratio of specific 

heats; finally, A(x) is the cross-sectional area of a quasi-one- 
dimensional duct. Written in the form (l), the Euler equations 
describe the propagation of pressure waves (here expressed in 

terms of Riemann variables, R^ and R 2 and entropy) along the 
characteristics having slopes X^ a nd X^, respectively, and the 
convection of entropy along streamlines. 

The interplay of pressure waves and entropy waves is an 

essential feature of unsteady flows. In a steady state, however, 
the third of (1 ) shows that the entropy must be constant wherever 
the flow is continuous and differentiable. Due to the presence 
of shocks in the flow, the entropy in a steady state is generally 
piecewise constant (note that contact surfaces cannot exist in a 
one-dimensional steady flow). If only the steady solution of a 
problem is of interest and we conceive of it as the asymptote of 
an unsteady process, the entropy terms can be dropped from (O, 
in order to find the steady solution in each region of continuous 
flow. This does not mean that the entropy does not influence the 
steady values of R 1 and R 2 * In fact, its influence is felt 
through the boundary conditions, as we will see further on. 

We will, consequently, simplify the Euler equations in two 
successive steps. In the first, we drop all entropy terms: 
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(5) 


It 


+ Mix + 


0 


2t 


+ X 2 R 2x + 


0 


The above system, which is correct only when the entropy is 
piecewise constant, will be used here solely for heuristic ar- 
guments, instead of the clumsier system (l). The following sys- 
tem, instead, is correct in all regions of continuous, steady 
flow, and it is the system whose numerical solution is sought: 


X 1 R 1x + d - 0 


X 2 R 2x + d * 0 


( 6 ) 


Let us assume that u is positive and consider subsonic flows 
only, for the time being. In subsonic flows, is negative and 
*2 is positive. The way of coding unsteady Euler solvers, known 
as the ^-scheme [l ], stems from the obvious fact, expressed by 
(5), that R 1 -waves propagate in the direction of decreasing x, 
whereas i^-waves propagate in the direction of increasing x. The 
two equations (5), thus, can be integrated separately at each 
computational node. Since the domains of dependence of R^ and R 2 
lie to the right and to the left of the point to be computed, 
respectively, the x-derivatives of R 1 and R 2 will have to be ap- 
proximated by finite differences between the point in question 
and points on its right or its left, respectively. 


Let t=kAt and x=nAx. To a first order of accuracy, (5) can 
be approximated by: 


R 

R 


k+1 
1 ,n 

k+1 

2,n 



0 tV?,..! - <,„> * d 


(7) 

( 8 ) 


Here, a = At /Ax. In so doing, we have laid the foundation for an 
explicit code. Alternatively, we can approximate (5) by: 


’f 1 = R* - Hf 1 ) + d A X ] 

1 ,n 1 ,n L 1 1 ,n+1 1 ,n J 


k+1 


,k+1 


’o +1 = Ro - R* +1 .) + d A X ] 

2,n 2,n L 2 2,n 2,n-1 J 


k+1 D k+1 


k+1 k+1 

Then R, and R~ are obtained as 
1 ,n ^ ,n 


(9) 

(10) 
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R, . + “ a 
1 > n 


[m h ?: 1 

TS 


R 1,n+1~ d Ax] 


1 +0)0 


( 11 ) 



JLjl 


+ ,lr d 

1 + u>aX„ 


Ax] 


(12) 


and an implicit code will be generated. In (11) and (12), a re- 
laxation factor, a) has been introduced, in the spirit of relaxa- 
tion techniques. 


Let us discuss the explicit approach first. In applying the 
explicit X-scheme to describe an unsteady flow, time must be ad- 
vanced by the same At at all computational nodes. The CFL cri- 
terion for stability [ 1 0 ] requires a to be the minimum of the 
values of 1/X^ at all points. The x advance in time is, thus, 
strongly limited. Conversely, in a computation which does not 
aim to provide a correct description of a time evolution, At can 
be let vary from point to point, provided that the discretized 
equations so obtained are compatible with the equations for 
steady flow, once the values of the physical parameters are lo- 
cally stabilized. 


Therefore, if we decide to integrate (7) and (8) separately 
over the entire region, and we begin by (8), we find that we can 
define a locally, as equal to I/X 2 , so that (8) becomes: 



R 


k 

2,n-1 


-d 


Ax/X£ 


(13) 


This formula can be used to integrate (8) without acting on (7), 
sweeping on the computational mesh from left to right. The pro- 
cedure is in the spirit of relaxation techniques, but it can 
still be interpreted as a description of the propagation of a 
right-running wave if the appearance of (13) as a particular case 
of (8) is taken into account. In the standard, explicit integra- 
tion of time-dependent equations, both right- and left-running 
waves are considered at every interval in order to update the 
flow correctly over all intervals but for a very short increment 
in time. Here, we attempt to describe the propagation of only 
one wave over the entire region without too much concern for the 
interaction of waves of the opposite family. Consequently, the 
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value of I ?2 at the right boundary point is, as far as the in- 
fluence of the right-running waves is concerned , updated over a 
time of the order of the length of the duct, divided by an 
average value of u+a. As we proceed in this phase of the in- 
tegration, we update u and a at every point, using the new value 
of R 2 and the old value of R^ , which is left unchanged. 

The procedure defined by (13) can be interpreted as an im- 
plicit calculation, as in (12), with an infinitely large value of 
w. The use of (13) is thus justified as optimal. In reaching 
the right boundary, we must apply the boundary condition. In 
problems of practical interest, the pressure is assigned at the 
exit of the duct. Insofar as there are no entropy sources, this 
is tantamount to assigning the speed of sound. Therefore, u can 
be evaluated from the second of (2) since R ^ has just been upda- 
ted. Then, R 1 follows from the first of (2). It is now time to 
integrate R^ , sweeping from right to left. Here we may not start 
from (7) and simplify it as we did with (8) since may get 
close to zero and the increment produced by d Ax/X^ may become 
too large and produce instabilities. We resort, instead, to 
(11), using a suitable value of to. Larger values of to speed up 
the convergence, but the theoretical gain rapidly tends to an 
asymptotic value. The advantage of the implicit formulation (ll) 
over the explicit one stems from the fact that (ll) is stable for 
large values of w even if X^ tends to zero. Again, as in the 
preceding sv iep, we will update u and a at every point, using the 
existing value of Rg and the new value of R 1 . 

In reaching the left boundary, we apply the boundary condi- 
tion by imposing the stagnation speed of sound, a Q . Since 

a 2 - a 2 + «u 2 (14) 


and R^ is accepted as computed, we obtain 


u 


-SR 


1 


+ 



S + 1 



(15) 


Then, a follows from (9) and R 2 from the second of (2). One in- 
tegration step is thus completed , and the computation resumes 
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with the right-running sweep. 

Note that (13) and ( 1 1 ) , regardless of the value of w, are 

compatible with the partial differential equations for the steady 

lc+1 

state, that is, when R!^ coincides with R* , to within the 
first order of accuracy. 


3. Second-order accurate scheme 


To achieve second-order accuracy in the asymptotic steady 
state solution, we approximate (6) by 


- <»> * <d n*1 + d „> 431 ‘ 0 

- " In -. 5 4 (d „-1 4 d »> 431 ‘ 0 


( 16 ) 


That such expressions are second-order accurate is easily seen. 
For example, the left-hand side of the first equation may be ex- 
panded, about point x=nAx, in the form: 


[2X + X’Ax + ^ X 1 -A x 2 ] [r'A x + ^ R"A X 2 ] + 
+ [2d + d' A x + -jr d’’Ax 2 ] A x 


(17) 


(all unessential indices having been dropped). By grouping terms 
with the same power of Ax and taking into account the first of 
(6), which we rewrite here: 


XR* + d = 0 


(18) 


and its x-derivative : 


X ' R ' + XR" + d’ = 0 (19) 

we see that (17) begins with terms of the third order in Ax. 
Consequently, (13) and (ll) are replaced by: 

‘ E 2,„-l - < d „-. 4d „> 4 */< X 2,„ 4 l 2,„-l> (?0 > 
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( 21 ) 


? k+1 
1 ,n 


R k + too 
1 » n 


+x. 


^ X 1,n'"1,n+1 

i +ua U , + x, 777 

1 1 ,n 1 ,n+1 1 


R^ 1 .. - (d +d .. )Ax] 
1 ,n+1 n n+1 J 


A FORTRAN code for the procedure outlined above is shown in 
Fig. 1. Most of the notation is obvious and requires no clarifi- 
cation; NC is the total number of nodes, B is the logarithmic in- 
crement of the area of the duct, A^/A.Ax), ANC is the prescribed 
speed of sound at the exit of the duct, Y is the stagnation speed 
of sound, in the nondimensionalization adopted here, GD,GE,GG,GC 
and GDH mean 5, 1+5, 1/5, (l+5)/5 and 5/2, respectively, and OM 
is wa. One may note the extreme speed of the calculation. Each 
double sweep requires only 12(NC-l) multiplications for interior 
points, and 15 multiplications for the boundary conditions, a to- 
tal of 12.NC+5 multiplications, as opposed to 58.NC+14 multipli- 
cations per step, as in a standard ^-scheme procedure, which is, 
per se, a fast-working code. Such an acceleration per step of 
almost a factor of 4 is compounded by a much higher speed of con- 
vergence, as we will show in what follows. 


£. First examples . Subsonic flows 

To illustrate the procedure, we present three calculations 
of subsonic flows in ducts. In all cases, the geometry of the 
duct is obtained by prescribing the steady-state Mach number dis- 
tribution. 

In the first case, we consider a convergent duct with a 
linear Mach number distribution between M*0.2 and M*0.8. We ap- 
ply the procedure, as coded in Fig. 1, using as initial condition 

1 /2 

a state of rest throughout the duct with a * Y , except for the 
exit endpoint, where a is set equal to the value belonging to the 
steady state solution. First, we check the accuracy. We repeat 
the computation five times, using 4, 8, 16, 52 and 64 intervals 
respectively. The mean quadratic error in u, after convergence 
is reached, is shown in Table I. The products of the second 
column by the square of the first column, reported in the third 
column, clearly show that the procedure is accurate to the second 
order. 
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SUBROUTINE FOR SUBSONIC FLOW ONLY 


N-1 

DNA«A(l )*U(1 )*B( 1 ) 

AL2A-U(l )+A(l ) 

DO 1 N»2,NC 
DN-A(N)*U(N)*B(N) 

AL2«U(N)+A(N) 

R (2 , N )«R (2 , N-1 )- (DN+DNA )/( AL2+AL2A ) 
U(N)«.5*(R(2,N)-R(1,N)) 

A(N)-GDH*(R(2 f N)+R(l ,N)) 

DNA-DN 

1 AL2A-AL2 

U (NC )=U (NC )+ (A(NC )-ANC )*GG 
A(NC)=ANC 

R(1 ,NC)-A(NC)*GG-U(NC) 

DNA=A(NC)*U(NC)*B(NC) 

AL1A«U(NC)-A(NC) 

DO 2 NN=2,NC 

N-NC+1-NN 

DN-A(N)*U(N)*B(N) 

AL1 -U (N )-A(N ) 

AL1 B=ABS (AL1 +AL1 A ) 

R(1 ,N)-(R(1 ,N)+0M*(AL1B*R(1 ,N+1 )-DN-DNA )/( 1 . +OM*AL1B) 
U (N )■=. 5*(R (2 , N)-R(l ,N)) 
a(n)-gdh*(r(2,n)+r(i ,N)) 

DNA-DN 

2 AL1A-AL1 

U(1 )=(-GD*R(l ,1 )+SQRT(GC*GAMMA-GD*R(2, 1 )**2))/GE 
A ( 1 )-SQRT (GAMMA -GD*U(1 )**2) 

R ( 1 , 1 )-GG*A(l )-U(l ) 

R(2,1 )=GG*A(1 )+U(l ) 

RETURN 

END 


Fig 1 . Subroutine for subsonic flow only 
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N 

E 

N 2 E 

WCT 

work 

E 

4 

1 .118E-3 

0.01 79 

1 

1190 

6.845E-5 

8 

2.81 8E-4 

0.01 79 

2 

697 

7.01 IE-5 

16 

7.038E-5 

0.0180 

4 

442 

7.051E-5 

32 

1 .659E-5 

0.01 70 

8 

544 

7.058E-5 

64 

2.830E-6 

0.01 16 

16 

1921 

7.082E-5 

N*=number of intervals 

32 

935 

7.034E-5 

E=mean 

quadratic error 

in u 

64 

1836 

7.085E-5 


TABLE I 



TABLE 

II 


Next, we test the speed of convergence. To this effect, we 
consider as a unit of work the work needed to process one mesh 
point in the two sweeps of each computational step. Therefore, 
the total work for one computation is the product of the tota] 
number of points (NC) by the total number of steps. We recompute 
the flow using 16 intervals and increasing values of mo (l, 2, 4, 
8, 16, 32 and 64). Table II shows the work needed to reach con- 
vergence. In the computer used in these tests (PDP 11/23)» we 
accept that convergence is reached when the mean square differen- 
ce in u between two successive steps is less than t=5x 10 For 
values of wo > 4, however, we had to increase t to 10" . In the 
same Table, we also report the mean quadratic error in u at con- 
vergence, to evidence that the results are not affected by the 
choice of It may be noted that wa=4 i s an optimum value. 
Increasing it not only hampers the convergence but creates oscil- 
lations in the results from step to step; a disturbing 
phenomenon, even if the oscillations are only of the order of 
10 ~ 7 . 


At this stage, we can try to compare the efficiency of the 
present calculation and of a standard time-dependent A-scheme 
technique. The comparison, however, offers some difficulties, 
since the converged accuracy of the standard technique which we 
used is worse than the converged accuracy of the present techni- 
que, when the same number of intervals is used. For example, for 
the same data and the same number of intervals, the standard 
technique converges in about 250 steps. Since the number of mul- 
tiplications needed at every point, that is the work per point, 
is about 4 times as large , we need thus a total equivalent work 
of the order of 17000. The economical gain (a factor of 40) is 
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LOG t 


remarkable, but it is grossly underestimated. In fact, the mean 
quadratic error in the standard calculation is about 7x10”^; with 
the present technique it is about 7*10 . To reach a similar ac- 
curacy with the standard technique ^t it3 present stage , we would 
need more points and, consequently, a much larger amount of work. 
For example, the work needed to reach an accuracy of 1.75*10“^ 
with the standard technique is equal to 80000 (almost 200 times 
as large as with the new technique), and the accuracy is still a 
quarter of an order of magnitude lower. We can tentatively con- 
clude that the new technique is more efficient than the standard 



Fig 2. Mean square residuals and errors - Convergent duct 


one by three orders of magnitude . 

Finally, we can reduce the total work even more by working 
with a multiple grid. Since, in this case, we do not need more 
than 16 points to achieve a very high accuracy, we will use first 
a set of 8 intervals and then a set of 16, with an automatic 
switch when convergence is reached in the coarser mesh. In going 
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from a coarser mesh to a finer one, we use second-order interpo- 
lations. The total work (in the case, wa= 4 ) equals 423. Mean 
square residuals and mean square errors are plotted in Fig. 2 for 
this case; work is on the abscissa. In Fig. 3 and 4 isobars are 
plotted for the present technique with x on the abscissa and work 
on the ordinate, and for the standard calculation with x on the 
abscissa and time on the ordinate. It is clear that the gain in 
convergence is achieved to the expense of an accurate description 



Fig 3. Isobars in the (x,work)-plane - Convergent duct 
Accelerated technique 


of the unsteady evolution. 

The second case is similar to the first but only apparently 
as simple. The duct is divergent, with a linear Mach number dis- 
tribution between 0.8 and 0.2. We use the same code and the same 
type of initial conditions as in the preceding case. Convergence 
is much harder to reach now, regardless of the computational 
technique. For conciseness, we present only results obtained 
with a multiple grid (8 and 16 intervals), and increasing values 
of uo. Our results are summarized in Table III. Mean square 
residuals and mean square errors in u are plotted in Fig. 5 for 
the case uo= 64; work is on the abscissa. We have also tried a 
run with three levels of mesh refinement (8, 16 and 32 intervals, 
successively), and one with four levels, up to 64 intervals. The 
value of T in all these cases was 5 X 10 . In the three-level 
run, a total work of 5497 was needed to reach E=1.06*10 - ^. In 
the four-level run, the computation stopped after a work of 6797, 
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Fig. 4. Isobars in the (x,t)-plane - Convergent duct 
Standard technique 


with E*5 • 88x 1 0“^ . The theoretical accuracy of 2.65x10“^ for a 
second-order accurate technique was not reached , but for this one 
needs a computer working on more digits and a lower value of t. 
As far as the standard technique is concerned, we have run a case 
with the same number of intervals and reached convergence after a 
work equal to 160000 but with an accuracy of only 6x10 . Using 
32 intervals, we reached an accuracy of 1.8*10 with a total 
work equal to 528000. 


(DO 

work 

E 

1 

14850 

5.48E-4 

2 

8107 

4.96E-4 

4 

5049 

4.96E-4 

8 

4265 

4.49E-4 

16 

3857 

4 . 49E-4 

32 

5311 

4.49E-4 

64 

3154 

4.49e-4 


TABLE III 




Again, the efficiency of the present technique is at least three 
orders of magnitude better than the efficiency of the standard 
technique. For the sake of comparing the acceleration in conver- 
gence during the initial phase of the calculation, in Fig. 6 and 
7 isobars are plotted for the present technique with x on the 
abscissa and work on the ordinate, and for a standard calculation 
with x on the abscissa and equivalent work (that is, number of 




Fig. 5. Mean square residuals and errors - Divergent duct 


points times number of steps times 4) on the ordinate. 

Finally, we consider a third case, in which a symmetric, 
parabolic Mach number distribution is assigned, with initial and 
final values of 0.2 and a maximum value of 0.9. The correspon- 
ding duct is now a symmetric nozzle with a throat. In this case, 
it took a work equal to 15406 to reach an accuracy of 3.44*10 , 
using a mesh of 16 intervals initially and then, for very few 
steps, a mesh of 32 intervals. The value of uor in this case was 
2. In Fig. 8 the exact distribution of the Mach number is shown 
by a solid line, and the computed values are denoted by 0's. 
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WORK 


10290 



Pig. 6. Isobars in the ( x,vork)-plane - Divergent duct 
Accelerated technique 



Fig. 7. Isobars in the (x,t)-plane - Divergent duct 
Standard technique 
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Fig. 8. Exact and computed M(x) for a symmetric nozzle 
with a throat 

5 _. Technique for transonic flows 

In Section 2, it was explained how the sweeping integration 
concept originated from the physical evidence of right-running 
I^-waves and left-running -waves. Consistency of numerics and 
physics in the handling both of interior points and of boundary 
conditions is, in our opinion, responsible for the high efficien- 
cy of the code, both in speed and accuracy. In this respect, it 
may be interesting to note that the technique described above is 
a hybrid between explicit and implicit techniques. It seems that 
the wave character of the unsteady flow has to be preserved in 
order to maintain a steady pattern which, per se, does not show 
an explicit propagation of waves; minute perturbations are easily 
carried away to the boundaries and properly adjusted. 

Such considerations should be kept in mind when attempting 
an extension of the technique to transonic flows. Here we 
describe a procedure for a case where a supersonic region is 
comprised between two subsonic regions. The interest of the case 
resides not only in its appearance in Laval nozzles with shocks, 
but in its close relation with patterns about supercritical air- 
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foils 


A few considerations are in order, before outlining a tech- 
nique . 

1 ) We assume that the supersonic region is bounded by a 
shock downstream. This is, indeed, the only case of practical 
interest. 

2) In a steady configuration, the entropy has only two 

values; one to the left of the shock (which can be assumed as 

zero) and one to the right of the shock. 

5) The R 2 -waves always propagate from left to right. 

4) The R 1 -waves propagate from right to left in subsonic 
regions and from left to right in the supersonic region. 

5) R 2 is discontinuous across the shock. Its jump, 

however, is very small. Therefore, the value of R^ behind the 

shock may be computed applying (20), but using only local infor- 
mation for d and and then incrementing it by its shock jump, 
as given below. 

6) The shock makes the set of values of R^ on its left and 

the set of values of R^ on its right totally independent of one 

another, a fact which is automatically reflected by the code, 
when it is written to satisfy property 4) above. 

7) The jump in R^ across the shock, AR^ f divided by the 

speed of sound in front of the shock, a ,, is a single-value 
function of the relative Mach number of the shock, defined by 

U s1 " W 

« * ( 22 ) 

si 

where u , is the velocity in front of the shock and W is the 

si 

shock velocity. Therefore, if a g1 and A Rj are known, M is also 
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known 


8) The sign of changes across the sonic point. Conse- 
quently, the updating of R 1 in the subsonic region at the left 
should start at the last subsonic point, without using informa- 
tion from the next supersonic point, and its updating in the su- 
personic region should start at the first supersonic point, 
without using information from the subsonic region (the informa- 
tion is transmitted by R^). 

The code for the transonic problem, shown in Pig. 9, is more 
complicated than the subsonic code shown in Fig. 1. The total 
number of operations, however, is not increased. A few comments 
will facilitate the reading of Fig. 9- The symbols, NSON and NS, 
denote the value of N at the first supersonic point and at the 
first point following the shock, respectively; DR2 is the jump of 
R2 across the shock. Loop 21 determines the location of the son- 
ic point at every step. The calculation of R^ starts at N=NC and 
proceeds backwards until the point denoted by NS has been compu- 
ted; then it jumps to the statement labeled 4, which resumes the 
calculation at the last subsonic point in the first region. The 
calculation proceeds backwards again until point N=1 is evalua- 
ted; then it jumps to the statement labeled 6, which resumes the 
calculation at the first supersonic point and proceeds forward 
until the last supersonic point has been computed. Finally, the 
shock is evaluated, the point behind it is corrected accordingly 
( in a separate subroutine) and the boundary condition at N=1 is 
applied as in the subsonic code. The exit boundary condition is 
imposed here by giving the logarithm of the exit pressure, PE 
( the pressure of the gas at rest being 1 ) and using the exit en- 
tropy, S(NC) as equal to the entropy produced by the shock. 


6_. Shock calculation 

Let us assume, for argument's sake, that we know in which 
interval the shock is located. Then, all points to the left of 
it are computed properly by our technique; indeed, R^ depends 
only on the left boundary condition and R^ is not influenced by 
the right boundary condition, due to the presence of the shock 
(recall property 6) above). To evaluate properly the points to 
the right of the shock, however, two elements are needed: the 
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SUBROUTINE FOR TRANSONIC FLOWS 


DNA«A(l )*U(1 )*B(1 ) 

AL2A-U(l )+A(l ) 

DO 1 N-2,NC 

dn=*a(n )*u (n )*b(n) 

DR2A-0. 

IF(N.NE.NS)GO TO 15 
DNA-DN 

AL2A*U(N)+A(N ) 

DR2A-DR2 
15 AL2-U (N )+A(N) 

DNB-DN+DNA 

AL2B-AL2+AL2A 

DNA-DN 

R(2,N)-R(2,N-1 )-DNB/ AL2B+DR2A 
U(N)-.5*(R(2,N)-R(1 * N) ) 
A(N)-GDH*(R(2,N)+rO ,N)) 

1 AL2A=AL2 
NSON1 -2 
NS0N2-NC 

IF(NSON.LE.NSON1+3)GO TO 22 
NSON1 “NSON-3 
NS0N2-NS0N+1 
22 DO 21 N-NSON1 , NS0N2 
21 IF(U(N).GE.A(N).AND.U(N-1 ) 
I.LT.A(N-I)) NSON-N 
ANC =SQRT (GAMMA*EXP (PE/GA +2 . 

1 *GD*S (NC ) ) ) 

U(NC)=U(NC)+(A(NC)-ANC)*GG 

A(NC)-ANC 

R ( 1 , nc)-a(nc )*gg-u(nc) 

DNA-A(NC)*U(NC)*B(NC) 

AL1 A-U (NC)-A(NC ) 

N-NC 
LITE--1 
3 H-N-1 

IF(N.EQ.NS-1 ) GO TO 4 


IF(N.LE.O) GO TO 6 
55 N1-N+1 

5 AL1»U(N)-A(N) 

DN-A(N)*U(N)*B(N) 

DNB-DNA+DN 
AL1B-AL1 +AL1A 
DNA-DN 

alib-abs(alib) 

R(1 ,N)«(R(l ,N)-0«*(-AL1B* 

1 R (1 ,N1 ) +DNB ) ) / ( 1 . +OM*AL1 B ) 
U(N)-.5*(R(2,N)-R(l t N)) 
A(N)-GDH*(R(2,N)+R(l ,N)) 
AL1A-AL1 

IF (LITE. EQ. 1 )GO TO 7 
GO TO 3 
4 N-NSON-1 

DNA-A(N)*U(N)*B(N) 

AL1A-U(N)-A(N) 

GO TO 55 

6 IF(NS. EQ. -1 0)G0 TO 13 
N-NSON-1 

DNA -A ( NSON )*U (NSON )*B(NSON ) 

AL1 A-U (NSON )-A(NSON ) 

LITE-1 

7 N-N+1 

IF(N.EQ.NS) GO TO 8 

N1-N-1 

GO TO 5 

8 CALL SHOCK 

13 U ( 1 )=(-GD*R(l ,1 )+SQRT(GC*GAMMA 
1 -GD*R (1,1 )**2 ) ) /GE 
A ( 1 )-SQRT (GAMMA-GD*U ( 1 )**2) 

R (2 , 1 )-GG*A( 1 )+U(l ) 

R(1 ,1 )-GG*A(l )-U(l ) 

RETURN 

END 


Fig. 9. Subroutine for transonic flows 
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value of at the point called NS, from which the values of I ?2 
at all the following points depend, and the value of S which, 
combined with the exit pressure, defines the exit speed of sound 
and, consequently, the proper exit boundary condition and R^ from 
the exit point to the point behind the shock. 

For a correct solution of the problem, thus, the shock must 

be fitted. Assuming, again, that the abscissa, x , of the shock 

s 

is known, R^ and R 2 can be extrapolated from upstream to the 
shock point in order to find the values of u and a at point A im- 
mediately in front of the shock (called u^ and a A , respectively), 
and R 1 can be extrapolated from downstream to find its value, 
R^ s , immediately behind the shock. Let M be the relative Mach 
number of the shock, as defined in (22), with u A and a A in lieu 
of u g ^ and a s1 , respectively; and let 

E - AR 1 /a A (23) 

Then, M and £ are related by the equation: 


= r(YM 2 -g)(l+6M 2 )] 1/2 + 6(M 2 -1 ) 1 

6 ( 1 +i)M « 


(24) 


which can be approximated, for 1_<M<4 t by the linear function: 


M = 1 + g] Z 

where g^ is defined by computing (25) at M=4: 

= 44(1 + 4) 

Sl 11«-4 + [(16 y-4)(1+16«)] 1 / 2 


(25) 


Therefore, M can be found with four correct digits by first using 
(25), then by computing a correction to £, using a secant method: 

- 2 £ - £' 

where £' is obtained from (24), and finally by using (25) again: 

M = 1 + g^, 
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Having so found the shock Mach number, three important quan- 
tities can be evaluated: the entropy behind the shock, 


S 


1 


[in 


vM 2 -<5 

1+6 


- V 


In 


(1+6)M ? I 
? 

1 +&n 


(26) 


the jump in across the shock: 


4R 2 ' »A » * 1 (27) 

and the shock velocity, W, given by (22). The new value of R 2 at 
N=NS is thus found by first getting R ? at point B immediately 
behind the shock, 


R 2B " R 2A 


+ Ar, 


(28) 


and then by interpolation: 


D R 2B + nR 2,NS+1 

2, NS = 1 + n 


(29) 


where 


h “ ( X JJg “ Xg)/^X (30) 

Having now both and R 2 at point N=NS, u and a can be recompu- 
ted at the same point. 


The shock velocity has no direct meaning, in a computation 
which aims to describe a steady state, where W is obviously zero. 
Nevertheless, it plays a very important role in moving the shock 
before convergence is reached. It is important to move the shock 
in order to locate it in the proper interval. To this effect, x g 
is updated at every step by adding W At to it, multiplied by a 
factor, a. The effect of such a factor is simply to accelerate 
the motion of the shock (a value of 3 seems to be convenient); it 
has obviously no effect on the steady state location of the 
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SHOCK SUBROUTINE 


EPS = (XS-X (NS-1 ))*DDX 
ETA=*1 . -EPS 

R 1 S =R ( 1 ,NS)+ETA*(R(l ,NS)-R(1 ,NS+1 )) 

AS1 =A(NS-1 )+EPS*(A(NS-1 )-A(NS-2 ) ) 

US1 *U (NS-1 )+EPS*(U (NS-1 )-U (NS-2 ) ) 

R2S 1 =GG*AS2+US2 
DR1 A=(R1S-R2S1 +US1 +US1 )/AS1 
EM*=1 . +GS1 *DR1 A 
SMIN=EH**2 

DR 1 B-2 . *DR 1 A- ( (SQRT ( (GAMMA*SMIN-GD )*( 1 . +GD*SMIN ) ) +GD* 

1 (SMIN-1 . ) )/(GE*EM)-1 . )*GG 
EM=1 . +GS1*DR1B 
SMIN=EM**2 
DUM=GE*EM 
DIM«=GD*SMIN 
DEM m GAMMA*SMIN-GD 
W=US1 -AS1 *EM 
XS-XS+W*DX*ALPHA 

SS2- ( ALOG (DEM/GE )-GAMMA*ALOG (GE*SMIN/( 1 . +DIM) ) )/(2 . *GD*GAMMA ) 
DR2«AS1 *(DR1 A +2. /(GE*EM)*(l . -SMIN) ) 

R2S2=R2S1 +DR2 

R (2 , NS )= (R2S2+ETA*R (2 , NS+1 ) )/(2 . +ETA ) 

U(NS)=.5*(R(2,NS)-R(1 ,NS)) 

A (NS )=GDH*(R (2 , NS )+R (l , NS ) ) 

IDUM=XS*NA 
NS=IDUM+2 
DO 1 N=NS,NC 
1 S(N)-SS2 

IF ( NS . EQ . NSO )RETURN 

if(ns.gt.nso)go TO 14 

R(1 ,NS)-2.*R(1 ,NS0)-R(1 ,NSO+1 ) 

R(2,NS)=2.*R(2,NS0)-R(2,NS0+1 ) 

GO TO 3 

14 R(1 ,NS0)=2.*R(l ,NSO-1 )-R(l ,NS0-2) 

R (2 , NSO )-2 . *R (2 , NSO-1 )-R (2 , NSO-2 ) 

3 NSO=NS 

U(NS)=.5*(R(2,NS)-R(1 ,NS)) 

a(ns)*gdh*(r(2,ns)+r(i ,ns)) 

RETURN 

END 


Fig. 10. Shock calculation 
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shock, since W vanishes. 

A code for the shock calculation is shown in Fig. 10. 
There, NSO is the value of NS at the preceding step, DDX*=1 /Ax and 
GS1“g^. The last part of the subroutine provides a smooth tran- 
sition when the shock point travels from one interval to the 
next, in either direction. 

We conclude by showing how the presence of a shock is detec- 
ted for the first time. So long as the flow is shockless, every 
interval is monitored by computing the pertinent value of £ and, 
the corresponding M, defined by (25). A shock is said to exist 
in the first interval for which M>1.1. Its abscissa is then de- 
fined as the middle point of the interval and the shock is track- 
ed, from that moment on, as described above. 


7_. Results for flows with shocks 

To test the code, we chose a nozzle so defined that, in the 
absence of shock, the Mach number would be linearly distributed 
between 0.5 and 2.5* Then we assume that a shock is present at 
x“0.63 (the total length of the nozzle being l). In such a way, 
the subsonic flow behind the shock can be evaluated, and the exit 
pressure to generate it is made known. The logarithm of the lat- 
ter is PE, mentioned in Section 5. In Table IV, we summarize 
results of runs made using 8, 16, 32 and 64 interval's. The total 
work, the shock location (x g ), the shock Mach number (m), and the 
mean quadratic error in u (E) are given. 


N 

work 

X 

g 

M 

E 

8 

900 

0.6244 

1.770 

3.09E-1 

16 

3400 

0.6256 

1.752 

5.32E-3 

32 

9180 

0.6281 

1.756 

3.11 e-4 

64 

16445 

0.6295 

1.757 

1 .36E-4 


TABLE IV 


Note that the exact location of the shock, as said above, is 
0.6300, and the exact shock Mach number is 1.7600. Another run 
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was made, taking 30 steps with 8 intervals, 50 steps with 16 in- 
tervals and then letting the number of steps for 32 and 64 inter- 

n 

vals be controlled by t*5x 10“ (for fewer intervals, due to 
sizeable oscillations of the shock, the mean quadratic residual 
never drops below 10”^). After a work equal to 16445, the compu- 
tation stops with the shock located at 0.6285, and a shock Mach 
number equal to 1.7571. The mean quadratic error in u is now 
1 .37x10” . In Figs. 11 and 12 we show the mean quadratic error 
and the mean quadratic residual for the latter case, as functions 
of the work. In Figs. 13 through 16 we show the exact and compu- 
ted Mach number distribution, entropy distribution, mass flow and 
stagnation speed of sound squared (as before, exact values are 
shown by solid lines and computed values by 0's). The only graph 
showing a flaw is the one for the mass flow. Errors generated at 
the shock are equivalent to a mass source. It must be noted, 
however, that the error in the mass flow is only 2/10 of 1$. 
Finally, Fig. 17 shows how fast the correct shock location is 
reached (x g is on the ordinate; work is on the abscissa). 

Once more, we wish to compare the speed of convergence of 
the present technique with the speed of convergence of a standard 
technique. The results are presented in Figs. 18 and 19; Fig. 18 
shows the isobars for the present run and Fig. 19 shows the iso- 
bars for a really time-dependent run. The work spent in reaching 
the top of Fig. 19 is equal to 156000 (one order of magnitude 
larger than with the present technique); but the mean quadratic 
error in u is 7.74*10“^ (almost two orders of magnitude larger). 

It should be noted that minor changes in the coding, which 
do not produce different truncation errors but merely different 
round-off errors, may affect the location of the shock and its 
Mach number, or the mass flow and the total temperature 
downstream of the 3hock by about 1/10 of \%. For example, if 
(27) is substituted by the equivalent expression: 

AR 2 = S T l ' +O ' M ( [ (yM 2 - 6)(1+6M 2 )] 1/2 - (1 +6 )M + 6 _ 6 M 2 ) (31) 

and the last run is repeated, the converged distributions of mass 
flow and total temperature are as shown in Figs. 20 and 21 . It 
is, indeed, hard to draw a conclusion from experiments of this 
kind. It is clear that, upstream of the shock, the accuracy can 
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LOGARITHMIC SCALE 



Fig. 11. Mean square error for flow with shock - Multiple grid 


be controlled as simply and efficiently as in the subsonic cases. 
One can see from the isobar plots (Fig. 18, for example) that 
minor variations at the shock upset the entire flow downstream of 
it. Any attempt to fix the problem, however, produces results 
which are exposed to the same degree of randomness as the calcu- 
lations mentioned above, because this is in the nature of things. 
In practice, we may conclude that, so long as an error less than 
1/10 of \% in the shock or behind it is irrelevant, all our 
results are identical, regardless of minor changes in the coding, 
changes in factors such as wo, a or B, and changes in the initial 
conditions. 

We have also conducted tests using different Mach number 
distributions (that is, different geometries) and different shock 
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Fig. 15. Exact and comput 
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locations, and we have found similar results. For example, with 
the same geometry as above, we obtained the results summarized in 
Table V. 


shock location shock Mach number 


theor. 

comp. 

theor. 

comp. 

0.6000 

0.5999 

1 .700 

1 .700 

0.6100 

0.6095 

1 .7200 

1 .712 

0.6200 

0.6218 

1 . 7400 

1.741 

0.6250 

0.6258, 

1 .7500 

1 .748 


TABLE V 



Fig. 18. Isobars in the (x, work) -plane - Nozzle with shock 
Accelerated technique 


B. Conclusions 

Ve have presented a technique for the numerical treatment of 
one-dimensional, steady, inviscid flows. The technique is ex- 
tremely fast and accurate. The total labor involved in the cal- 
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TIME 



Pig. 19. Isobars in the (x,t)-plane - Nozzle with shock 

Standard technique 



Pig. 20. Exact and computed aass flow for flow with shock 


- 30 - 


WORK 




1*2 


MttZ 
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Inin m i m il i n n 111 1 11 11 1 1 1111 1 1 111111 i n ri 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ri 

Fig. 21. Exact and computed a^(x) for flow with shock 

culation (and, consequently, the total running time) is at least 
three orders of magnitude smaller than for a standard, time- 
dependent technique with the same degree of accuracy. Shocks are 
fitted with great accuracy also. It is important to note that 
our estimates of accuracy are not made by comparing numerical 
results of runs made with different number of intervals, but by 
comparing the computed results to exact solutions. Therefore, 
our "errors" are real errors, and should not be confused with 
what is call "error", at times, by other authors. The latter is 
instead what we called "residual", that is, the mean square dif- 
ference between values at one step and values at the preceding 
step. 


One final word should be spent on the subject of resolution. 
We have seen, indeed, that good results in flows with shocks are 
obtained when a relatively large number of nodes is used (in our 
examples above, we have used a maximum of 64 intervals), but 
reasonable results are also obtained with fewer intervals, as Ta- 
ble IV shows. Some deterioration on local values of sensitive 
parameters, such as mass flow and total temperature, may be ex- 


- 31 




HS4? , RUH *> 1 H YMIN»YHAX* 0.00000E 00 0.20000E 01 



NS47 » RUN 101 HASS YHIN,YMAX= 0.U800E 01 0.12500E 01 



Fig. 22. Exact and computed M(x) and mass flow for flow with shock 

8 intervals 
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NS47 , RUN 102 


YHINi YHAX= 0.00000E 00 0.20000E 01 



HS47 < RUN 102 NASS YN1N.YHAX- 0.11500E01 012000E 01 



Fig» 23- Exact and conputed H(x) and Bass flow for flow with shock 

16 intervals 
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NS47 , RUN 183 NASS YHIN.VHAX* 0.11500E01 0.12000E 01 



Fig. 24. Exact and coaputad «(x) and mass flow for flow with shock 

32 intervals 
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pected if the number of computational nodes is very small. It is 
important to develop a feeling for what is more likely to 
deteriorate, in view of extensions to two-dimensional problems, 
where the number of mesh points is often much smaller than in 
one-dimensional cases. To this effect, we present plots of M and 
PuA as functions of x, in Figs. 22 through 24, having used 8, 16 
and 32 intervals respectively. The Mach number is perfect in all 
cases; the mass flow behind the shock, instead, deteriorates when 
few interval are used. 
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